FoRTE baseline

Abstract

Forests in the US Upper Midwest have been a net C sink for a century, driven by their fast-growing, early successional species. However, the composition of these regrown forests has started to shift toward mid-successional species. Combined with changing climate and disturbance, this means that future biogeochemical and ecological dynamics in these forests are highly uncertain. Demographic vegetation models provide a tool for analyzing this uncertainty and its implications. Here, we evaluate the relative importance of model “structural” (which processes are represented, and how) and “parameter” (parameters used in those equations) uncertainty, asking: (1) Which processes are most important to accurately model community succession and C cycling in Upper Midwest forests? (2) What are the relative contributions of parameter vs. structural uncertainty? We ran ensembles of the Ecosystem Demography model v2.2 (sampling over parameter uncertainty) from 1900 to 2000 with different representations of key processes important to competition for light. We then compared the magnitude of structural (across representations) and parameter (within each representation) uncertainty, and assessed which submodel-parameter combinations reproduced observed C fluxes and succession. Our simulations agreed well with observed net primary productivity but systematically underestimated leaf area index. Overall, parameter uncertainty contributed more than structural uncertainty to total simulation uncertainty, but the amount of parameter uncertainty varied with model structure. We conclude that: (1) targeted data collection to constrain vegetation model parameters should be a higher priority than model development; (2) parameter sensitivity analyses are not necessarily generalizable across model structures; and (3) robust model validation requires scrutiny of multiple different variables, including those related to vegetation productivity, structure, and composition.

Introduction

Forests are a defining feature of the US North. This region is the most heavily forested region in the US, and these forests are an essential provider of employment and recreation opportunities and other ecosystem services for the region’s people, and of habitat for the region’s wildlife (Shividenko et al. 2005, Shifley et al. 2012). In addition, the forests of the Upper Midwest in particular have been a net carbon sink for most of the past century, as fast-growing early successional species established and thrived following near-complete deforestation around the end of the 19th Century (Birdsey et al. 2006, Williams et al. 2012). However, the fate of these forests is highly uncertain. For one, many Upper Midwest forests are undergoing a shift in community composition from early-successional to mid-successional species, and the ecological and biogeochemical consequences of this transition are not well-established (Gough et al. 2016). In addition, these forests face a multitude of direct and indirect pressures from both biotic and abiotic sources (Shifley and Moser 2016), including longer and more severe droughts (Gustafson et al. 2016, Liénard et al. 2016, Swanston et al. 2017), non-native insects and pathogens (Lovett et al. 2016), and more. Collectively, these processes have culminated in forests that are different from their pre-settlement climax counterparts (Wolter and White 2002, Stearns and Likens 2002, Thompson et al. 2013), and our ability to make predictions about such non-analog environments based on the past is limited.

More reliable predictions can in theory be obtained by using dynamic vegetation models that explicitly represent processes involved in forest growth and mortality. Vegetation models fall broadly along the following spectrum of complexity, which is a direct consequence of trade-offs between realism and process fidelity on one hand and tractability, computational demand, and data requirements on the other (Hawkes 2000). On one side are relatively simple “big leaf” models—such as PNET, SiBCASA, and Biome-BGC—in which the vegetation at a particular location consists essentially of a single large “plant” whose characteristics are the (weighted) average of all the vegetation at that site. A common application of these models was simulating the land surface boundary condition for atmospheric general circulation models, though similar models have been (and continue to be) applied successfully to simulate biogeochemical processes in many ecosystem science contexts. On the other end of the spectrum are “individual-based” models (a.k.a. “gap models”), which explicitly simulate multiple indivduals competing for resources at a single site (Shugart et al. 2015). Examples of such models are LANDIS (Mladenoff 2004, Scheller et al. 2007) and UVAFME. Because these models can explicitly represent inter-specific differences in plant productivity, resource allocation, and stress tolerance (among others), as well as the competition that emerges out of these differences, they may be better able to represent changes in ecosystem-scale processes, especially in no-analog conditions (Purves and Pacala 2008). Between these two extremes are models that use approximations attempt to capture the emergent biogeochemical and ecological outcomes of interactions between individual plants without the need for explicitly simulating each individual (Purves and Pacala 2008). One specific example of this approach is the Ecosystem Demography model (ED2; Moorcroft et al. 2001, Medvigy and Moorcroft 2011, Longo et al. 2019b), which is the focus of this study.

Regardless of vegetation model complexity, their projections are inherently uncertain. This uncertainty comes from many different sources, which can broadly be classified into the following categories: Driver uncertainty refers to uncertainty in data about processes not represented by the model (e.g. weather and climate for vegetation models). Initial condition uncertainty arises from the fact that models have to start somewhere, and the exact conditions at the place and time simulations begin are frequently unknown. Process or structural uncertainty arises because vegetation models necessarily only represent a subset of all processes involved in plant biology, and that the process that are included can usually be represented in multiple different ways (i.e. which processes are included, and how they are represented). Finally, parameter uncertainty arises because of natural variability and imperfect calibration of the above process representations. The relative importance of these sources of uncertainty in ecological forecasts is a key question because it is directly related to future research priorities. For example, in the atmospheric science community, the insight that meteorological forecasts were most sensitive to initial condition uncertainty (Lorenz REF) led to a multi-decadal research agenda aimed at constraining the initial state with improved observations, which has directly contributed to a steady and persistent improvement in meteorological forecasts over the last several decades.

Many ecosystem modeling studies have looked at these categories of uncertainty independently. Perhaps the most work has been done on parameter uncertainty in ecosystem models. Dietze et al. (2014) investigated which parameters were most important to ED2 predictions of C sequestration across a diverse range of sites in North America. The found that the parameters contributing the most to overall predictive uncertainty were those related to growth respiration, mortality, stomatal conductance, and water uptake, and that parameter uncertainty varied with biome. More recently, Raczka et al. (2018) evaluated parameter uncertainty over 100 year time scale at the Willow Creek Ameriflux site, and found that the most important parameters to model uncertainty were quantum efficiency of photosynthesis, leaf respiration, and soil-plant water transfer. Another recent study assessed parameter uncertainty in ED2 as part of its overall goal of parameterizing the model for sagebrush (Pandit et al. 2018). They found that the most important parameters were specific leaf area, the maximum rate of carbon fixation (V_cmax), slope of stomatal conductance, the fine-root to leaf carbon ratio, and the fine root turnover rate. Another study looked specifically at contributions to uncertainty from parameters related to canopy radiative transfer, and found that parameters related to both canopy structure and leaf optical properties had a large impact on predictions of net primary productivity (Viskari et al. 2019). {Driver uncertainty – a few examples?} Other studies have investigated structural uncertainty.. Numerous vegetation model intercomparison studies have demonstrated that different models produce significantly different projections of overall land carbon sequestration (Friedlingstein et al. 2006, 2014), response to CO₂ fertilization (Zaehle et al. 2014, Walker et al. 2015, Medlyn et al. 2015), and soil C sequestration (Sulman et al. 2018), among others. In particular, Lovenduski and Bonan (2017) found that model structural uncertainty was the primary source of uncertainty in Earth system model projections of the land C sink. This structural uncertainty has been attributed to differences in model representations of key processes, including canopy radiative transfer (Fisher et al. 2017), soil biogeochemistry (Wieder et al. 2017, Sulman et al. 2018), stomatal conductance (Franks et al. 2018), and photosynthesis (Rogers et al. 2016). However, the extent to which specific processes contribute to model uncertainty is difficult to evaluate from model intercomparisons because different models are different from each other in too many different ways, meaning that structural effects are confounded with other uncertainty sources. Moreover, comparative analysis of contributions of different types of uncertainty to model projections are rare in the ecosystem modeling literature.

This study focuses on the interaction between parametric and structural uncertainty in the ED2 model. Our study is organized around the following guiding questions: (1) Which processes related to light utilization are most important to accurately modeling community succession and C cycle dynamics in temperate forests of the Upper Midwest? (2) What is the cost of considering these processes, in terms of additional parametric uncertainty? (3) What are the relative contributions of parametric uncertainty (data limitation) vs. structural uncertainty (theoretical limitation?)? To answer these questions, we ran ED2 ensemble simulationes with a factorial combination of submodels related to radiative transfer formulation (two-stream vs. multiple scatter), horizontal competition (finite canopy radius vs. complete shading), and trait plasticity (whether or not SLA and Vcmax vary with light level) We then used sensitivity and variance decomposition analyses to evaluate the contribution of parameter uncertainty for each model configuration, exploring the ecological implications of the results.

Methods

Site description

We performed this study at the University of Michigan Biological Station (UMBS; Ameriflux site US-UMd, 45.5625\(^{\circ}\), -84.6975\(^{\circ}\)), located in Northern Lower Michigan, USA. The area surrounding the research station is 87% well-drained upland forest and 13% wetland (Bergen and Dronova 2007); the focus of this study is on the former. The landscape geography of the UMBS upland forest is 20.4% moraine, 37.8% high outwash plain, 31.3% low outwash plain, 5.7% lake-margin terrace, 3.6% ice-contact, and 1.2% lowland glacial lake (Bergen and Dronova 2007). Most of the UMBS upland forest canopy is dominated by temperate deciduous early-successional species, most importantly Populus grandidentata (bigtooth aspen) and, to a lesser extent, Betula papyrifera (paper birch), with Acer saccharum (sugar maple), Acer rubrum (red maple), Fagus grandifolia (American beech), Tilia americana (basswood), Betula alleghaniensis (yellow birch), Fraxinus americana (white ash), Tsuga canadensis (eastern hemlock), Quercus rubra (northern red oak), Pinus strobus (white pine), and Pinus resinosa (red pine) existing in various fractions in the understory (and, in patches, in the canopy). This composition is a legacy of the site’s disturbance history: The site was intensively logged in the late 1800s and early 1900s, and experienced regularly recurring fires until the mid-1920s, at which point a regime of active fire suppression started that has persisted to the present day. As a result, the average stand age in 2013 was 95 years. The majority of the forest that is aspen-dominated is undergoing succession to “northern hardwood” (maple, beech, basswood, birch, ash, hemlock), “upland conifer” (red and white pine), or “northern red-oak” ecotypes (Bergen and Dronova 2007).

Forest stands at UMBS were previously exposed to experimental disturbance as part of the Forest Accelerated Succession Experiment (FASET). {Canonical FoRTE citation (Curtis and Gough 2018)}

Ecosystem Demography Model (version 2.2)

For this study, we used the Ecosystem Demography Model, version 2.2 (ED-2.2). A full description of the default configuration of the model is provided by Longo et al. (2019b). Briefly, ED-2.2 solves the energy, water, and carbon cycles separately for each of multiple “cohorts” of trees of similar composition, size, and age sharing a micro-environment and disturbance history (Moorcroft et al. 2001, Medvigy et al. 2009, Longo et al. 2019b). Besides the cohort-based approximation of gap dynamics, a distinctive feature of ED-2.2 is that the equations for energy, water, and carbon fluxes are defined in terms of total energy, water, and carbon, which ensures excellent conservation of mass in long-running (multidecadal) simulations (Longo et al. 2019b). Various versions of the ED model have been validated in boreal (Medvigy and Moorcroft 2011), temperate (Medvigy et al. 2009, Dietze et al. 2014, Raczka et al. 2018), and tropical biomes (Moorcroft et al. 2001, Longo et al. 2019a), and have been applied in, among others, paleoecological studies (Rollinson et al. 2017), free-air CO₂ enrichment studies (De Kauwe et al. 2013, Miller et al. 2015), estimation of potential carbon stocks (Hurtt et al. 2004), and analyses of regional vegetation-climate feedbacks (Swann et al. 2015).

The official ED2 source code is available at https://github.com/EDmodel/ED2, and the exact version used in this study (which includes minor revisions to accommodate the requirements of this study) can be obtained at https://github.com/ashiklom/ED2/tree/b048950971e91699b78fc566df133caf977fda32.

Submodels

In this study, we ran a factorial combination of the following ED-2.2 configurations: (1) two-stream vs. multiple-scatter canopy radiative transfer models; (2) infinite vs. finite crown area (a.k.a. complete vs. partial shading); and (3) static vs. light-plastic traits.

Both canopy radiative transfer models in ED-2.2 resolve the full vertical radiation profile within a patch as a function of canopy structure (leaf and wood area indices, crown area, leaf angle distribution) and incident solar radiation, following the definitions in CLM 4.5 (Oleson et al. 2013). Both models also use identical definitions for the optical properties of a single canopy layer as a function of leaf and wood optical properties, leaf angle distribution, canopy clumping, solar zenith angle, and leaf, wood, and crown area indices. Both models represent direct (a.k.a. “beam”) radiation as an exponentially decaying process, and solve for the diffuse (a.k.a. “hemispherical”, “isotropic”) radiation at each layer using a linear matrix equation. Where the models differ is in the terms of this matrix equation. The two-stream model (default) uses the two-stream approximation originally derived for atmospheric radiative transfer (Meador and Weaver 1980) and later adapted to vegetation canopies (Dickinson 1983, Sellers 1985). Meanwhile, the multiple-scatter model was derived from first principles by Zhao and Qualls (2005) specifically for vegetation canopies, specifically to address known biases and limitations of the two-stream approach (e.g. Wang 2003). Both of these models are described in detail in the Supplementary Information.

The crown area submodel in ED-2.2 determines the nature of competition for light between cohorts. In the default configuration (“infinite crown area”, or “complete shading”), the leaf area of a cohort is distributed across the entire horizontal area of a patch. This means that taller cohorts always receive more incoming radiation than shorter cohorts, even when the height difference is small. This has been shown to excessively suppress competition from sub-dominant individuals and result in unrealistically homogeneous canopies (Fisher et al. 2015). In the alternate configuration (“finite crown area”, or “partial shading”), canopies take up only a fraction of the available horizontal area, meaning that multiple cohorts of similar height can receive the same level of light. The horizontal area of crowns is determined by allometric equations from Dietze and Clark (2008).

The third submodel we evaluated was trait plasticity. In the default configuration, all cohorts of a given plant functional type will have the same parameters, regardless of environmental conditions. This ignores the globally-documented intraspecific trait variability as a function of light level (Niinemets 2010, Keenan and Niinemets 2016). In the alternate configuration, as light level decreases (trees become more shaded), specific leaf area increases and V_cmax decreases, following empirical relationships from the tropics (Lloyd et al. 2010).

These sub-models are summarized in Table 1.

Table 1: ED-2.2 submodel descriptions

Name	Description	Color
Crown model
closed (default)	Cohort crowns take up entire patch area. Competition for light based only on height.	Light
finite	Cohort crown area is proportional to DBH according to PFT-specific allometry.	Dark
Radiative transfer model
two-stream (default)	Two-stream approximation (Sellers 1985, Oleson et al. 2013)	Primary (red, blue)
multi-scatter	Multiple-scatter approximation, following (Zhao and Qualls 2005)	Secondary (green, orange)
Trait plasticity
static (default)	SLA and Vc,max are constant	Cool (blue, green)
plastic	SLA increases, and Vc,max decreases, with light level	Warm (red, orange)

Parameter descriptions

The full list of parameters in ED-2.2 is large, with over 100 parameters per plant functional type. A full sensitivity and uncertainty analysis across all of these parameters is outside the scope of this study. Instead, we selected a subset of parameters guided by (but expanding upon) previous ED2 sensitivity studies (Dietze et al. 2014, Raczka et al. 2018, Viskari et al. 2019).

Three parameters are related to leaf-level physiology: Following the enzyme-kinetic model of Farquhar et al. (1980), the rate of photosynthesis is the minimum of light-limited and enzyme-limited reactions. The former are controlled by the quantum efficiency parameter—maximum efficiency with which absorbed photosynthetically active radiation is converted to CO₂. The latter are controlled by V_cmax, the maximum rate of carbon fixation by Rubisco. The water demand of photosynthesis is controlled by the stomatal slope, the sensitivity of stomatal conductance of CO₂ as a function of CO₂ concentration and humidity at the leaf surface (Ball et al. 1987, Leuning 1995). Three more parameters correspond to the respiration rates of leaves, roots, and “growth maintenance”.

Two more parameters control carbon allocation: One is the ratio of fine root to leaf biomass (fineroot2leaf, or q), and another is the ratio of “storage” carbon allocated to reproduction (r_fract).

Three parameters control various aspects of adult tree mortality. The overall adult mortality rate in ED-2.2 (\(M\)) is the sum of density-independent mortality from aging (\(M_I\), plants year^-1), density-dependent mortality from C starvation (\(M_D\)), and mortality from cold/frost (\(M_F\)) (ED-2.2 technically also includes additional term for fire mortality, but we did not include fire in our simulations):

\[ M = M_I + M_D + M_F \]

Density-independent mortality from aging (\(M_I\)) is a prescribed, PFT-specific parameter (mort3). Density-dependent mortality from C starvation is calculated as a function of a cohort’s C balance limitation:

\[ M_D = \frac{y_1}{1 + \exp\left[y_2 \left( \frac{C_k}{C_k^*} \right)\right]} \]

where \(C_k\) is the 12 month running mean C balance of the cohort, \(C_k^*\) is the running mean ideal C balance if there was no light or water limitation, and \(y_1\) (mort1; plants year^-1) and \(y_2\) (mort2; unitless) are PFT-specific parameters. Seedling mortality rate is prescribed as its own PFT-specific parameter.

Several parameters are related to canopy structure and radiative transfer. Specific leaf area (SLA) is used to convert leaf biomass to leaf area index, which in turn is used in a variety of calculations related to canopy radiative transfer and micrometeorology. Canopy clumping factor describes how evenly leaf area is distributed in horizontal space (1 being perfectly evenly; 0 being a “black hole” where all leaves are concentrated in a single point); and leaf orientation factor describes the average distribution of leaf angles (-1 being perfectly vertical, 1 being perfectly horizontal, and 0 being random). Four parameters control leaf optical properties, namely the fractions of light reflected or transmitted in visible and near-infrared wavelengths (leaf_(reflect|trans)_(vis|nir)). Details of how these parameters influence canopy radiative transfer are described in the Supplementary Information.

Other parameters:

• Fraction of litter carbon that enters the labile (fast) soil carbon pool (f_labile)
• Minimum plant height for reproduction (minimum_height)
• Soil-plant water conductance, or “water availability factor” (water_conductance)
• C:N ratio in fine roots (c2n_fineroot) and leaves (c2n_leaf)
• Leaf turnover rate (evergreen-only; leaf_turnover_rate)
• Proportion of dispersal that is “global” (nonlocal_dispersal)

The full list of parameters is shown in Table 2.

Table 2: Parameter descriptions. For additional information, see the ED2 model wiki.

ED Name	Description	unit_markdown	Raczka	Dietze
c2n_fineroot	C:N ratio in fine roots	unitless (mass ratio)	NA	NA
c2n_leaf	C:N ratio in leaves	unitless (mass ratio)	NA	NA
clumping_factor	Canopy clumping factor	NA	NA	NA
f_labile	Fraction of litter that goes to labile (fast) C pool	unitless (0-1)	Labile carbon	NA
fineroot2leaf (q)	Ratio of fine root to leaf biomass	unitless (mass ratio)	Root/Leaf carbon	Leaf:Root
growth_resp_factor	Fraction of daily C gain lost to growth respiration	unitless (0-1)	Growth respiration	Growth Resp
leaf_reflect_nir	Leaf reflectance in NIR range (700-2500 nm)	unitless (0-1)	NA	NA
leaf_reflect_vis	Leaf reflectance in visible range (400-700 nm)	unitless (0-1)	NA	NA
leaf_respiration_rate (Rd0)	Ratio of leaf respiration to Vcmax	???	Leaf respiration*	Leaf Resp*
leaf_trans_nir	Leaf transmittance in NIR range (700-2500 nm)	unitless (0-1)	NA	NA
leaf_trans_vis	Leaf transmittance in visible range (400-700 nm)	unitless (0-1)	NA	NA
leaf_turnover_rate	Temperature dependent rate of leaf loss (conifer only)	year-1	NA	NA
minimum_height	Minimum height for plant reproduction	m	Minimum height	NA
mort1	Time-scale at which low-carbon balance plants die	years-1	Carbon balance mortality	Mortality
mort2	C balance ratio at which mortality rapidly increases	unitless (mass ratio)	NA	NA
mort3	Density-independent (background) mortality rate	year-1	Background mortality	NA
nonlocal_dispersal	Proportion of dispersal that is global	unitless (0-1)	NA	NA
orient_factor	Leaf angle orientation distribution	unitless (-1-1)	NA	NA
quantum_efficiency	Farquhar model parameter (TODO)	mol CO2 (mol photons)-1	Quantum efficiency	Quantum Eff.
r_fract	Fraction of C storage to seed reproduction	unitless	Recruitment carbon	Reproduction?
root_respiration_rate (_factor)	Root respiration rate at 15 °C	μmol CO2 (kg fineroot)-1	Root respiration	NA
root_turnover_rate	Temperature dependent rate of fine root loss	year-1	Root turnover	Root turnover
seedling_mortality	Proportion of seed that dies and goes to litter pool	unitless (0-1)	NA	NA
SLA	Specific leaf area	m2 kg-1 C	Specific leaf area	SLA
stomatal_slope	Slope between A and stomatal conductance (Leuning)	unitless	Stomatal sensitivity	Stomatal Slope
Vcmax (Vm0)	Maximum rate of CO2 carboxylation at 15 °C	μmol m-2 s-1	Vcmax	Vcmax
water_conductance	Water availability factor	m-2 a-1 (kg C root)-1	Soil-plant water conductance	Water Cond

Plant functional types and parameterization

By default, ED-2.2 supports 17 different plant functional types, which divide plant species according to photosynthetic pathway (C3 vs. C4), growth form (grass vs. tree), leaf phenology habit (deciduous vs. evergreen), biome (e.g. temperate vs. tropical), and successional status (e.g. early, mid, late). However, we limited our simulations to the four plant functional types that have any appreciable presence at UMBS: Early, mid, and late temperate deciduous trees and pines. The species comprising these plant functional types are shown in Table 3.

Table 3: Plant functional type-species mappings

Plant functional type	Species	Color
Early hardwood	Betula papyrifera	Violet
	Populus grandidentata
	Populus tremuloides
Mid hardwood	Quercus rubra	Blue
	Acer rubrum
	Acer pensylvaticum
Late hardwood	Acer saccharum	Green
	Fagus grandifolia
Pine	Pinus strobus	Yellow

For each plant functional type, we generated a distribution of parameter values via the Predictive Ecosystem Analyzer (PEcAn) trait-meta analysis (LeBauer et al. 2013, see also Dietze et al. 2014, Raczka et al. 2018). Prior distributions for this meta-analysis were largely adapted from previous ED2 parameter uncertainty studies (Dietze et al. 2014, Raczka et al. 2018, Viskari et al. 2019), and are shown in detail in Supplementary Information X. Species trait data for this meta-analysis came from existing records in the BETY database (www.betydb.org, LeBauer et al. 2017), as well as from publicly available records in the TRY database (www.try-db.org, Kattge et al. 2011) and, specifically for leaf optical properties, from Shiklomanov (2019 {dissertation, chapter 3}). The resulting parameter distributions are shown in Figure 1. All parameters not described in this section were set to their ED-2.2 PFT-specific default values.

Figure 1: Input parameter distributions from PEcAn trait meta-analysis.

Modeling workflow

For each factorial combination of ED-2.2 configurations (described above), we ran 100 ensemble members from 1901 to 2000. Each ensemble member was initialized from a “near-bare ground” condition: An equal number of seedlings of each plant functional type (see previous section) at the minimum resolvable size. Driving meteorological data was 6-hourly CRU-NCEP combined with an annual atmospheric CO₂ record from Law-Dome ice core (Etheridge et al. 1998) and Mauna Loa observatory (Thoning et al. 1989). Soil texture was set to 92% sand, 7% silt, and 1% clay, per site-level observations in (Gough et al. 2010). The initial soil moisture profile was set to the average soil moisture profile reported in the UMBS Ameriflux ancillary data (https://ameriflux.lbl.gov/sites/siteinfo/US-UMd).

Driver uncertainty is outside the scope of this study, and initial condition uncertainty is minimized by our experimental design, which is conditioned on a specific initial condition (near bare ground).

Analysis of results

• Overall summary:
  • Overall parameter uncertainty – Time-series plots of total annual GPP and NPP, and maximum LAI and AGB, averaged across ensemble members
  • Impact of parameter uncertainty on ecological trajectories – plot of max annual LAI by PFT
  • Both compared against observations of recently observed LAI and NPP at UMBS (e.g. FASET, FoRTE)
• Uncertainty analysis:
  • Parameter vs. structural uncertainty – Variance in mean annual GPP from 19XX to 19YY across ensembles within model structure, and compare against
  • Breakdown of parameter uncertainty – PEcAn sensitivity and uncertainty analysis (LeBauer et al. 2013, Dietze et al. 2014)
• Details of structural uncertainty:
  • Vertical profiles of radiation, traits?

All analyses for this work were performed using R 3.6 (see “Colophon” for full details). Code and supporting data for reproducing this analysis are publicly available on the Open Science Framework (OSF) at https://osf.io/dznuf/.

Results

Summary

Figure 2: ED2 plot-level predictions of gross (GPP) and net (NPP) primary productivity, aboveground biomass (AGB), total leaf area index (LAI), and Shannon diversity index (of plant functional types) by model configuration. The solid line is the ensemble mean. The dashed line is the 80% confidence interval. Black dot with error bars is the mean and min/max value from Hardiman et al. (2013).

Model projections of gross and net primary productivity, aboveground biomass, leaf area index, and functional diversity varied with both parameter value and model structure (Figure 2). One model configuration – finite canopy radius, two-stream radiative transfer, and plastic traits – failed to grow at all, regardless of parameter combinations, and a similar configuration with static traits showed much less growth than any other model configurations. The remaining configurations all exhibited a similar trajectory of rapid growth in the first 15-20 years followed by a gradual decline to a lower, stable value around 1975 (with individual ensemble members sometimes undergoing sudden collapse in specific years). The ensemble mean ecosystem state of most model confugrations had net primary productivity very close to that observed at UMBS, but generally much lower leaf area index. However, the ensemble spread around the mean was large – across all configurations (except the finite canopy - two-stream radiative transfer cases), the spread in net primary productivity was from zero to several times higher than the observed. Ensemble variance in leaf area index was lower on average but, in some case (e.g. finite canopy, multiple-scatter radiative transfer, static traits) was over an order of magnitude. The predicted Shannon diversity index of plant functional types was, on average, higher in simulations when the canopy architecture was finite rather than closed.

Figure 3: Post-1975 growing season averaged NPP vs. LAI, with observation (in black) and linear regression, by model type. Letter labels are selected runs with the same parameters, and match the rows in Figure 4.

In most simulations, the leaf area index was too low for a given level of net primary productivity, but the relationship between these two variables varied with model configuration (Figure 3). The crown model (closed vs. finite) had the strongest effect: Compared to the simulations with “closed” crowns, simulations with finite canopy radius generally had lower net primary productivity for the same leaf area index. Only a few individual simulations matched both the observed leaf area index and net primary productivity, and all of them had the finite canopy radius configuration.

Figure 4: PFT-level predictions of leaf area index (LAI) by model configuration and parameterization. Each row of plots is a simulation with the same set of parameters across model configurations. Letters correspond to the labels in Figure 3. Note that these are predictions of total leaf area index summed across all cohorts of the same PFT, meaning that higher LAI in one PFT does not necessarily indicate dominance, size, etc.

Model structure and parameters interacted to determine the community assembly and dynamics of each simulation, with consequences for overall ecosystem productivity that varied by model structure (Figure 4). Consider simulations A and C, the former showing competition of early hardwood with pine, and the latter completely dominated by pine: With a closed canopy and static traits, both simulations produced similar overall LAI and NPP values (Figure 3), but moving to a finite canopy dramatically increased the total LAI in C without much affecting A, while enabling trait plasticity dramatically increased the NPP of A without much affecting C (Figure 3). This example also shows the impacts of different configurations on predicted community dynamics – namely, in this specific case, enabling trait plasticity helped pine out-perform early hardwood, but only under a closed (rather than finite) canopy. Finite canopy radius similarly benefits early hardwood in its competition against late hardwood in simulations B and, to a lesser extent, D. Simulations B and D also show that the sensitivity of a given parameter combination to model structure is itself dependent on the parameters – simulation B produced similar ecological dynamics and aggregate state variables regardless of model configuration (except the extreme cases; Figure 3), while simulation D had near-zero NPP and LAI in two cases (closed, two-stream, static and finite, multi-scatter, plastic), comparable NPP and LAI to B with closed canopy and trait plasticity, and higher NPP and LAI than D with a finite canopy, multi-scatter radiative transfer, and static traits (Figure 3).

Uncertainty analysis

Figure 5: Comparison of parametric uncertainty within ensembles (colored bars, colored by model type) and “structural” uncertainty (variance in ensemble means; black bar) by output variable. Output variables are expressed as growing season averages for all years after 1910.

Except for the configurations with persistently low productivity, parameter uncertainty (across ensembles within each model configuration) was much larger than structural uncertainty (variance in ensemble means across configurations; Figure 4). The relative importance of parameter uncertainty varied with both model configuration and the variable in question. Allowing SLA and V_cmax to vary with light level (warm colors, as opposed to static traits – cool colors) slightly reduced variability in primary productivity, slightly increased variability in Shannon diversity, and had mixed effects on aboveground biomass and total leaf area. Enabling the finite canopy radius (dark colors; compared to closed canopy – light colors) dramatically decreased variability in primary productivity and Shannon diversity and increased variability in leaf area, and had mixed effects on aboveground biomass. reduced variance in aboveground biomass, diversity, and primary productivity. Finally, using a multiple-scatter (secondary colors) rather than a two-stream (primary colors) radiative transfer scheme increased variance in leaf area index, aboveground biomass, and diversity, but did not have a consistent effect on variance in productivity productivity.

Figure 6: PEcAn-like parameter sensitivity and uncertainty analysis, by model type, for 1975-2000 averaged growing season net primary productivity (left) and leaf area index (right). Elasticity (top) is the normalized sensitivity of the model to a fixed change in the parameter. Partial variance (bottom) describes the overall contribution of the parameter to model predictive uncertainty based on the combination of parameter uncertainty and model sensitivity. Each panel shows the top 10 parameter values, in terms of the absolute value of the corresponding quantity. Input parameter distributions are shown in Figure 1. Note that x-axis scales vary across panels.

Model sensitivity to specific parameters (“elasticity” sensu LeBauer et al. 2013) varied depending on the model configuration and the target output variable (Figure 5). The parameter with the single highest elasticity to NPP and LAI across all model configurations was specific leaf area (SLA), though the PFT whose SLA was most important varied with model configuration. Other frequently (but not always) important parameters were the leaf C:N ratio (c2n_leaf), ratio of fine root to leaf biomass (fineroot2leaf), leaf and root respiration rates, and leaf optical properties (leaf_reflect_vis/nir).

The parameters that contributed most to overall model predictive uncertainty (“partial variance”) differed from parameters to which the model was most sensitive (Figure 5), largely driven by the fact that many of the parameters with the highest sensitivity (e.g. SLA, V_cmax) were also well constrained by data (Figure 1). Leaf optical properties of at least one PFT were among the top 10 most important parameters across all configurations, and in multiple cases, were the most important. Other often important parameters were plant-soil water conductance (water_conductance), parameters related to adult mortality (mort1,2,3), and growth respiration. Again, patterns in partial variance were largely idiosyncratic across model configurations.

Discussion

• Parameter uncertainty is the dominant driver of model predictive uncertainty.

• …but the nature of this parameter uncertainty depends on model configuration.

• Ecological implications of structural uncertainty

  • Predictions about ecological succession are very sensitive to the understory light environment.
    • Need a good representation of canopy structure (finite canopy radius)
    • How you model canopy radiative transfer matters! Small differences in the model can give significant results.
  • But, trait plasticity can provide an alternative mechanism for success in the understory.
    • Important modeling implication for all the emerging research on intra-specific trait variability.
    • This combination of discrete functional types with plastic traits is closer to reality than purely trait-based (“eco-evolutionary?”) models or completely discretized models with static traits on the other.
    • Important caveat: Current implementation is based entirely on tropical data, which has more extreme canopy shading than temperate forests (REF), so tropical species may have evolved stronger trait plasticity than temperate species (REF).
    • Perhaps note “ecological equifinality” here: there are different ways to arrive at similar outcomes in the model; this is also true in nature?

• Implications of uncertainty analysis for future data collection and model development

  • Parameter uncertainty is still a huge source of overall model predictive uncertainty. Different parameters can give qualitatively different ecological outcomes.
  • Looking beyond aggregate variables (e.g. net primary productivity, total LAI) is essential to model evaluation and calibration. It’s possible to generate a relatively?
  • (Expected) Subsurface parameters (e.g. leaf-root ratios) have high sensitivity, but are poorly constrained.
  • Trait covariance may alleviate some data limitation (Shiklomanov et al. 2018 [preprint])

• Comparisons against other studies

  • Parameter sensitivity
    • (Pandit et al. 2018) – Most important parameters for shrubs were SLA, Vcmax, stomatal slope, root:leaf ratio, and root turonver rate (slightly less important: leaf turnover rate, growth respiration, cuticular conductance; much less important: water conductance; basically irrelevant: storage turnover, leaf width)
    • Compare against Dietze et al. (2014) and Raczka et al. (2018)

• Future directions: Predicting disturbance response

  • FoRTE project – Forecasting a disturbance.
    • Compare ecological resilience/memory between model and reality

• Other relevant references:

  • Gough et al. (2019) – Canopy structural complexity is the strongest predictor of productivity across NEON sites. Possible future direction?

Moreover, these models’ predictions of species composition and ecosystem structure are relevant outside of carbon cycle science, particularly for forest management and wildlife conservation.

Conclusions

Acknowledgements

This project funded by NSF grant. Cyberinfrastructure provided by Pacific Northwest National Laboratory (PNNL). Data from University of Michigan Biological Station (UMBS). Data from TRY (TODO: Specific TRY statement).

pagebreak

References

Ball, J. T., I. E. Woodrow, and J. A. Berry. 1987. A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions. Progress in Photosynthesis Research:221–224.

Bazzaz, F. A. 1979. The physiological ecology of plant succession. Annual Review of Ecology and Systematics 10:351–371.

Bergen, K. M., and I. Dronova. 2007. Observing succession on aspen-dominated landscapes using a remote sensing-ecosystem approach. Landscape Ecology 22:1395–1410.

Birdsey, R., K. Pregitzer, and A. Lucier. 2006. Forest carbon management in the united states. Journal of Environment Quality 35:1461.

Curtis, P. S., and C. M. Gough. 2018. Forest aging, disturbance and the carbon cycle. New Phytologist 219:1188–1193.

De Kauwe, M. G., B. E. Medlyn, S. Zaehle, A. P. Walker, M. C. Dietze, T. Hickler, A. K. Jain, Y. Luo, W. J. Parton, I. C. Prentice, and et al. 2013. Forest water use and water use efficiency at elevated CO2: A model-data intercomparison at two contrasting temperate forest FACE sites. Global Change Biology 19:1759–1779.

Dickinson, R. E. 1983. Land surface processes and climate-surface albedos and energy balance. Theory of Climate, Proceedings of a Symposium Commemorating the Two-Hundredth Anniversary of the Academy of Sciences of Lisbon:305–353.

Dietze, M. C., and J. S. Clark. 2008. Changing the gap dynamics paradigm: Vegetative regeneration control on forest response to disturbance. Ecological Monographs 78:331–347.

Dietze, M. C., S. P. Serbin, C. Davidson, A. R. Desai, X. Feng, R. Kelly, R. Kooper, D. LeBauer, J. Mantooth, K. McHenry, and et al. 2014. A quantitative assessment of a terrestrial biosphere model’s data needs across North American biomes. Journal of Geophysical Research: Biogeosciences 119:286–300.

Etheridge, D. M., L. P. Steele, R. L. Langefelds, R. J. Francey, J.-M. Marnola, and V. I. Morgan. 1998. Historical CO2 records from the Law Dome DE08, DE08-2, and DSS ice cores. In trends: A compendium of data on global change. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, TN, USA.

Farquhar, G. D., S. von Caemmerer, and J. A. Berry. 1980. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149:78–90.

Fisher, R. A., C. D. Koven, W. R. L. Anderegg, B. O. Christoffersen, M. C. Dietze, C. E. Farrior, J. A. Holm, G. C. Hurtt, R. G. Knox, P. J. Lawrence, J. W. Lichstein, M. Longo, A. M. Matheny, D. Medvigy, H. C. Muller-Landau, T. L. Powell, S. P. Serbin, H. Sato, J. K. Shuman, B. Smith, A. T. Trugman, T. Viskari, H. Verbeeck, E. Weng, C. Xu, X. Xu, T. Zhang, and P. R. Moorcroft. 2017. Vegetation demographics in earth system models: A review of progress and priorities. Global Change Biology 24:35–54.

Fisher, R. A., S. Muszala, M. Verteinstein, P. Lawrence, C. Xu, N. G. McDowell, R. G. Knox, C. Koven, J. Holm, B. M. Rogers, and et al. 2015. Taking off the training wheels: The properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED). Geoscientific Model Development 8:3593–3619.

Franks, P. J., G. B. Bonan, J. A. Berry, D. L. Lombardozzi, N. M. Holbrook, N. Herold, and K. W. Oleson. 2018. Comparing optimal and empirical stomatal conductance models for application in Earth system models. Global Change Biology 24:5708–5723.

Friedlingstein, P., P. Cox, R. Betts, L. Bopp, W. von Bloh, V. Brovkin, P. Cadule, S. Doney, M. Eby, I. Fung, and et al. 2006. Climate-carbon cycle feedback analysis: Results from the C4MIP model intercomparison. Journal of Climate 19:3337–3353.

Friedlingstein, P., M. Meinshausen, V. K. Arora, C. D. Jones, A. Anav, S. K. Liddicoat, and R. Knutti. 2014. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. Journal of Climate 27:511–526.

Friedman, S. K., and P. B. Reich. 2005. Regional legacies of logging: Departure from presettlement forest conditions in northern Minnesota. Ecological Applications 15:726–744.

Gates, F. C. 1930. Aspen association in northern Lower Michigan. Botanical Gazette 90:233–259.

Gough, C. M., J. W. Atkins, R. T. Fahey, and B. S. Hardiman. 2019. High rates of primary production in structurally complex forests. Ecology.

Gough, C. M., P. S. Curtis, B. S. Hardiman, C. M. Scheuermann, and B. Bond-Lamberty. 2016. Disturbance, complexity, and succession of net ecosystem production in North America’s temperate deciduous forests. Ecosphere 7:e01375.

Gough, C. M., B. S. Hardiman, L. E. Nave, G. Bohrer, K. D. Maurer, C. S. Vogel, K. J. Nadelhoffer, and P. S. Curtis. 2013. Sustained carbon uptake and storage following moderate disturbance in a great lakes forest. Ecological Applications 23:1202–1215.

Gough, C. M., C. S. Vogel, B. Hardiman, and P. S. Curtis. 2010. Wood net primary production resilience in an unmanaged forest transitioning from early to middle succession. Forest Ecology and Management 260:36–41.

Gough, C. M., C. S. Vogel, H. P. Schmid, H.-B. Su, and P. S. Curtis. 2008. Multi-year convergence of biometric and meteorological estimates of forest carbon storage. Agricultural and Forest Meteorology 148:158–170.

Gustafson, E. J., A. M. G. De Bruijn, B. R. Miranda, and B. R. Sturtevant. 2016. Implications of mechanistic modeling of drought effects on growth and competition in forest landscape models. Ecosphere 7:e01253.

Hawkes, C. 2000. Woody plant mortality algorithms: Description, problems and progress. Ecological Modelling 126:225–248.

Hurtt, G. C., R. Dubayah, J. Drake, P. R. Moorcroft, S. W. Pacala, J. B. Blair, and M. G. Fearon. 2004. Beyond potential vegetation: Combining lidar data and a height-structured model for carbon studies. Ecological Applications 14:873–883.

Kattge, J., S. Díaz, S. Lavorel, I. C. Prentice, P. Leadley, G. Bönisch, E. Garnier, M. Westoby, P. B. Reich, I. J. Wright, J. H. C. Cornelissen, C. Violle, S. P. Harrison, P. M. Van Bodegom, M. Reichstein, B. J. Enquist, N. A. Soudzilovskaia, D. D. Ackerly, M. Anand, O. Atkin, M. Bahn, T. R. Baker, D. Baldocchi, R. Bekker, C. C. Blanco, B. Blonder, W. J. Bond, R. Bradstock, D. E. Bunker, F. Casanoves, J. Cavender-Bares, J. Q. Chambers, F. S. Chapin III, J. Chave, D. Coomes, W. K. Cornwell, J. M. Craine, B. H. Dobrin, L. Duarte, W. Durka, J. Elser, G. Esser, M. Estiarte, W. F. Fagan, J. Fang, F. Fernández-Méndez, A. Fidelis, B. Finegan, O. Flores, H. Ford, D. Frank, G. T. Freschet, N. M. Fyllas, R. V. Gallagher, W. A. Green, A. G. Gutierrez, T. Hickler, S. I. Higgins, J. G. Hodgson, A. Jalili, S. Jansen, C. A. Joly, A. J. Kerkhoff, D. Kirkup, K. Kitajima, M. Kleyer, S. Klotz, J. M. H. Knops, K. Kramer, I. Kühn, H. Kurokawa, D. Laughlin, T. D. Lee, M. Leishman, F. Lens, T. Lenz, S. L. Lewis, J. Lloyd, J. Llusià, F. Louault, S. Ma, M. D. Mahecha, P. Manning, T. Massad, B. E. Medlyn, J. Messier, A. T. Moles, S. C. Müller, K. Nadrowski, S. Naeem, Ü. Niinemets, S. Nöllert, A. Nüske, R. Ogaya, J. Oleksyn, V. G. Onipchenko, Y. Onoda, J. Ordoñez, G. Overbeck, W. A. Ozinga, S. Patiño, S. Paula, J. G. Pausas, J. Peñuelas, O. L. Phillips, V. Pillar, H. Poorter, L. Poorter, P. Poschlod, A. Prinzing, R. Proulx, A. Rammig, S. Reinsch, B. Reu, L. Sack, B. Salgado-Negret, J. Sardans, S. Shiodera, B. Shipley, A. Siefert, E. Sosinski, J.-F. Soussana, E. Swaine, N. Swenson, K. Thompson, P. Thornton, M. Waldram, E. Weiher, M. White, S. White, S. J. Wright, B. Yguel, S. Zaehle, A. E. Zanne, and C. Wirth. 2011. TRY – a global database of plant traits. Global Change Biology 17:2905–2935.

Keenan, T. F., and Ü. Niinemets. 2016. Global leaf trait estimates biased due to plasticity in the shade. Nature Plants 3.

Kilburn, P. D. 1960. Effects of logging and fire on xerophytic forests in northern Michigan. Bulletin of the Torrey Botanical Club 87:402.

LeBauer, D., R. Kooper, P. Mulrooney, S. Rohde, D. Wang, S. P. Long, and M. C. Dietze. 2017. BETYdb: A yield, trait, and ecosystem service database applied to second-generation bioenergy feedstock production. GCB Bioenergy 10:61–71.

LeBauer, D. S., D. Wang, K. T. Richter, C. C. Davidson, and M. C. Dietze. 2013. Facilitating feedbacks between field measurements and ecosystem models. Ecological Monographs 83:133–154.

Leuning, R. 1995. A critical appraisal of a combined stomatal-photosynthesis model for C3 plants. Plant, Cell and Environment 18:339–355.

Liénard, J., J. Harrison, and N. Strigul. 2016. US forest response to projected climate-related stress: Atoleranceperspective. Global Change Biology 22:2875–2886.

Lloyd, J., S. Patiño, R. Q. Paiva, G. B. Nardoto, C. A. Quesada, A. J. B. Santos, T. R. Baker, W. A. Brand, I. Hilke, H. Gielmann, and et al. 2010. Optimisation of photosynthetic carbon gain and within-canopy gradients of associated foliar traits for Amazon forest trees. Biogeosciences 7:1833–1859.

Longo, M., R. G. Knox, N. M. Levine, A. L. S. Swann, D. M. Medvigy, M. C. Dietze, Y. Kim, K. Zhang, D. Bonal, B. Burban, and et al. 2019a. The biophysics, ecology, and biogeochemistry of functionally diverse, vertically- and horizontally-heterogeneous ecosystems: The Ecosystem Demography model, version 2.2 - part 2: Model evaluation. Geoscientific Model Development Discussions:1–34.

Longo, M., R. G. Knox, D. M. Medvigy, N. M. Levine, M. C. Dietze, Y. Kim, A. L. S. Swann, K. Zhang, C. R. Rollinson, R. L. Bras, and et al. 2019b. The biophysics, ecology, and biogeochemistry of functionally diverse, vertically- and horizontally-heterogeneous ecosystems: The Ecosystem Demography model, version 2.2 - part 1: Model description. Geoscientific Model Development Discussions:1–53.

Lovenduski, N. S., and G. B. Bonan. 2017. Reducing uncertainty in projections of terrestrial carbon uptake. Environmental Research Letters 12:044020.

Lovett, G. M., M. Weiss, A. M. Liebhold, T. P. Holmes, B. Leung, K. F. Lambert, D. A. Orwig, F. T. Campbell, J. Rosenthal, D. G. McCullough, and et al. 2016. Nonnative forest insects and pathogens in the United States: Impacts and policy options. Ecological Applications 26:1437–1455.

Meador, W. E., and W. R. Weaver. 1980. Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement. Journal of the Atmospheric Sciences 37:630–643.

Medlyn, B. E., S. Zaehle, M. G. De Kauwe, A. P. Walker, M. C. Dietze, P. J. Hanson, T. Hickler, A. K. Jain, Y. Luo, W. Parton, and et al. 2015. Using ecosystem experiments to improve vegetation models. Nature Climate Change 5:528–534.

Medvigy, D., and P. R. Moorcroft. 2011. Predicting ecosystem dynamics at regional scales: An evaluation of a terrestrial biosphere model for the forests of northeastern North America. Philosophical Transactions of the Royal Society B: Biological Sciences 367:222–235.

Medvigy, D., S. C. Wofsy, J. W. Munger, D. Y. Hollinger, and P. R. Moorcroft. 2009. Mechanistic scaling of ecosystem function and dynamics in space and time: Ecosystem demography model version 2. Journal of Geophysical Research 114.

Miller, A. D., M. C. Dietze, E. H. DeLucia, and K. J. Anderson-Teixeira. 2015. Alteration of forest succession and carbon cycling under elevated CO2. Global Change Biology 22:351–363.

Mladenoff, D. J. 2004. LANDIS and forest landscape models. Ecological Modelling 180:7–19.

Moorcroft, P. R., G. C. Hurtt, and S. W. Pacala. 2001. A method for scaling vegetation dynamics: The ecosystem demography model (ed). Ecological Monographs 71:557–586.

Niinemets, Ü. 2010. A review of light interception in plant stands from leaf to canopy in different plant functional types and in species with varying shade tolerance. Ecological Research 25:693–714.

Oleson, K. W., D. M. Lawrence, G. B. Bonan, B. Drewniak, M. Huang, C. D. Koven, S. Levis, F. Li, W. J. Riley, Z. M. Subin, S. C. Swenson, P. Thornton, A. Bozbiyik, R. Fisher, C. L. Heald, E. Kluzek, J.-F. Lamarque, P. J. Lawrence, R. Leung, W. Lipscom, S. Muszala, D. M. Ricciuto, W. Sacks, Y. Sun, J. Tang, and Y. Zong-Liang. 2013. Technical description of version 4.5 of the Community Land Model (CLM). NCAR Earth System Laboratory Climate; Global Dynamics Division.

Pandit, K., H. Dashti, N. F. Glenn, A. N. Flores, K. C. Maguire, D. J. Shinneman, G. N. Flerchinger, and A. W. Fellows. 2018. Optimizing shrub parameters to estimate gross primary production of the sagebrush ecosystem using the Ecosystem Demography (EDv2.2) model. Geoscientific Model Development Discussions:1–23.

Purves, D., and S. Pacala. 2008. Predictive models of forest dynamics. Science 320:1452–1453.

Raczka, B., M. C. Dietze, S. P. Serbin, and K. J. Davis. 2018. What limits predictive certainty of long‐term carbon uptake? Journal of Geophysical Research: Biogeosciences 123:3570–3588.

Rogers, A., B. E. Medlyn, J. S. Dukes, G. Bonan, S. von Caemmerer, M. C. Dietze, J. Kattge, A. D. B. Leakey, L. M. Mercado, Ü. Niinemets, and et al. 2016. A roadmap for improving the representation of photosynthesis in earth system models. New Phytologist 213:22–42.

Rollinson, C. R., Y. Liu, A. Raiho, D. J. P. Moore, J. McLachlan, D. A. Bishop, A. Dye, J. H. Matthes, A. Hessl, T. Hickler, and et al. 2017. Emergent climate and CO2 sensitivities of net primary productivity in ecosystem models do not agree with empirical data in temperate forests of eastern North America. Global Change Biology 23:2755–2767.

Scheller, R. M., J. B. Domingo, B. R. Sturtevant, J. S. Williams, A. Rudy, E. J. Gustafson, and D. J. Mladenoff. 2007. Design, development, and application of LANDIS-II, a spatial landscape simulation model with flexible temporal and spatial resolution. Ecological Modelling 201:409–419.

Sellers, P. J. 1985. Canopy reflectance, photosynthesis and transpiration. International Journal of Remote Sensing 6:1335–1372.

Shifley, S. R., F. X. Aguilar, N. Song, S. I. Stewart, D. J. Nowak, D. D. Gormanson, W. K. Moser, S. Wormstead, and E. J. Greenfield. 2012. Forests of the Northern United States. U.S. Department of Agriculture, Forest Service, Northern Research Station.

Shifley, S. R., and W. K. Moser. 2016. Future forests of the northern United States. U.S. Department of Agriculture, Forest Service, Northern Research Station.

Shividenko, A., C. V. Barber, R. Persson, P. Gonzalez, R. Hassan, P. Lakyda, I. McCallum, S. Nilsson, J. Pulhin, B. van Rosenburg, and B. Scholes. 2005. Forest and woodland ecosystems. in R. Hassan, R. Scholes, and N. Ash, editors. Millennium ecosystem assessment: Ecosystems and human well-being: Current state and trends, volume 1. Island Press.

Shugart, H. H., G. P. Asner, R. Fischer, A. Huth, N. Knapp, T. Le Toan, and J. K. Shuman. 2015. Computer and remote-sensing infrastructure to enhance large-scale testing of individual-based forest models. Frontiers in Ecology and the Environment 13:503–511.

Stearns, F., and G. E. Likens. 2002. One hundred years of recovery of a pine forest in northern Wisconsin. The American Midland Naturalist 148:2–19.

Sulman, B. N., J. A. M. Moore, R. Abramoff, C. Averill, S. Kivlin, K. Georgiou, B. Sridhar, M. D. Hartman, G. Wang, W. R. Wieder, and et al. 2018. Multiple models and experiments underscore large uncertainty in soil carbon dynamics. Biogeochemistry 141:109–123.

Swann, A. L. S., M. Longo, R. G. Knox, E. Lee, and P. R. Moorcroft. 2015. Future deforestation in the Amazon and consequences for South American climate. Agricultural and Forest Meteorology 214-215:12–24.

Swanston, C., L. A. Brandt, M. K. Janowiak, S. D. Handler, P. Butler-Leopold, L. Iverson, F. R. Thompson III, T. A. Ontl, and P. D. Shannon. 2017. Vulnerability of forests of the Midwest and Northeast United States to climate change. Climatic Change 146:103–116.

Thompson, J. R., D. N. Carpenter, C. V. Cogbill, and D. R. Foster. 2013. Four centuries of change in northeastern United States forests. PLoS ONE 8:e72540.

Thoning, K. W., P. P. Tans, and W. D. Komhyr. 1989. Atmospheric carbon dioxide at Mauna Loa Observatory: 2. Analysis of the NOAA GMCC data, 1974-1985. Journal of Geophysical Research: Atmospheres 94:8549–8565.

Viskari, T., A. Shiklomanov, M. C. Dietze, and S. P. Serbin. 2019. The influence of canopy radiation parameter uncertainty on model projections of terrestrial carbon and energy cycling. PLOS ONE 14:e0216512.

Walker, A. P., S. Zaehle, B. E. Medlyn, M. G. De Kauwe, S. Asao, T. Hickler, W. Parton, D. M. Ricciuto, Y.-P. Wang, and D. W. et al. 2015. Predicting long-term carbon sequestration in response to CO2 enrichment: How and why do current ecosystem models differ? Global Biogeochemical Cycles 29:476–495.

Wang, Y. P. 2003. A comparison of three different canopy radiation models commonly used in plant modelling. Functional Plant Biology 30:143.

Wieder, W. R., M. D. Hartman, B. N. Sulman, Y.-P. Wang, C. D. Koven, and G. B. Bonan. 2017. Carbon cycle confidence and uncertainty: Exploring variation among soil biogeochemical models. Global Change Biology 24:1563–1579.

Williams, C. A., G. J. Collatz, J. Masek, and S. N. Goward. 2012. Carbon consequences of forest disturbance and recovery across the conterminous United States. Global Biogeochemical Cycles 26:n/a–n/a.

Wolter, P. T., and M. A. White. 2002. Recent forest cover type transitions and landscape structural changes in northeast Minnesota, USA. Landscape Ecology 17:133–155.

Zaehle, S., B. E. Medlyn, M. G. De Kauwe, A. P. Walker, M. C. Dietze, T. Hickler, Y. Luo, Y.-P. Wang, B. El-Masri, P. Thornton, and et al. 2014. Evaluation of 11 terrestrial carbon-nitrogen cycle models against observations from two temperate free-air CO2 enrichment studies. New Phytologist 202:803–822.

Zhao, W., and R. J. Qualls. 2005. A multiple-layer canopy scattering model to simulate shortwave radiation distribution within a homogeneous plant canopy. Water Resources Research 41.

Appendix 1: Submodel descriptions

Radiative transfer: Definitions

• \(\omega\) – Clumping factor
• \(LAI\) – Leaf area index
• \(TAI\) – Total area index (leaf and woody area)
• \(LAI_e\) / \(TAI_e\) – “Effective” area indices (accounting for clumping)

Two-stream RTM

Multiple scatter RTM

Effective LAI is just true LAI x clumping.

\[ LAI_e = \omega LAI \]

Total area index (TAI) is leaf (LAI) + wood (WAI) area.

\[ TAI_e = LAI_e + WAI \]

Local exposed total area index (\(TAI_l\)) is the total area index divided by crown area index (\(CAI\)).

\[ TAI_{l} = \frac{TAI_e}{CAI} \]

Leaf (\(w_l\)) and wood (\(w_w\)) weights:

\[ w_l = \frac{LAI_e}{TAI_e} \] \[ w_w = 1 - w_l \]

Projected area (\(a_{proj}\)), based on coefficients \(\phi_1\) and \(\phi_2\).

\[ a_{proj} = \phi_1 + \phi_2 \mu \]

Optical depth (\(\lambda\)) is projected area over optical depth (\(\mu\))

\[ \lambda = \frac{a_{proj}}{\mu} \]

Transmittance for direct radiation of a canopy layer (\(\tau_{beam}\)) is:

\[ \tau_{beam} = (1 - CAI) + CAI exp\left( -\lambda TAI_l \right) \]

Transmittance for diffuse radiation of a canopy layer (\(\tau_{diff}\)):

\[ ext_1 = \phi_1 TAI_l \] \[ ext_2 = \phi_2 TAI_l \] \[ \tau_{diff} = (1 - CAI) - CAI exp\left( -ext_1 - ext_2 \right) ext_1^2 e^{exp_1} E_i(-exp_1) + (ext_1 - 1)\]

where \(E_i\) is the exponential integral.

Single-scattering albedo (\(w_0\)):

\[ w_0 = \frac{1}{2} \frac{a_{proj}}{\phi_2 \mu + a_{proj}} \frac{1 - \phi_1 \mu}{\phi_2 \mu a_{proj}} \log \left( 1 + \frac{\phi_2 \mu + a_{proj}}{\phi_1 \mu} \right) \]

Beam backscatter (\(\gamma\)):

\[ \gamma = w_0 \frac{1 + \bar{\mu} \lambda}{\bar{\mu} \lambda} \]

Old outline

Introduction

• Why we need models?
  • Forests are important, for C sequestration, other ecosystem services
  • Policy decisions related to forests depend on future projections
  • We need models to make future projections
• What kinds of models are out there?
  • Traditional: “Big leaf” models
    • Needed to be simple to run on old computers
    • Don’t need as much data
    • …but they fail to capture important processes, which matter not only for accurate predictions about C cycle, but also because they are processes we care about (habitat, forestry, tourism, etc.)!
  • Individual-based models
    • Much more detailed, but harder to run, and need much more data.
    • Trade-off between complexity and accuracy.
• Cohort-based models like ED2 as a middle ground?
  • Core ED concept (Moorcroft et al. 2001)
  • ED2 developments – (Medvigy et al. 2009, Medvigy and Moorcroft 2011)
  • Most recent version – (Longo et al. 2019b)
  • Validation/applications of the ED model:
    • (Rollinson et al. 2017) – ED2 paleoclimate
    • (Miller et al. 2015) – ED2 at Duke FACE
    • (De Kauwe et al. 2013) – ED2 at Duke and Oak Ridge FACE
    • (Longo et al. 2019a) – Model evaluation in Tropics
• How much complexity do we need to accurately model ecological dynamics?
  • Lots of existing work on parametric uncertainty
    • (Dietze et al. 2014)
    • (Raczka et al. 2018)
    • (Pandit et al. 2018)
  • Less on structural uncertainty – i.e. not just what coefficients do we use, but what equations, and more generally, what processes do we even consider?
• What are the processes that drive succession of ecological communities?
  • Competition for light (Bazzaz 1979)
  • Others? Leaning towards omitting because factorial combination becomes too large.
    • Competition for water?
    • Competition for nutrients? (ED’s N model isn’t very good. Alternatives are being developed, but not in ED mainline yet.)
• Brief review of community/ecosystem ecology of the Upper Midwest
  • UMBS is a useful site for testing these theories because of its natural history.
  • Around 1900, mature forest dominated by white (Pinus strobus) and red pine (Pinus resinosa), hemlock (Tsuga canadensis), and sugar maple (Acer saccharum).
  • But, repeat harvest and burning between 1880 and 1920 (Kilburn 1960) led to replacement by early successional aspen (Populus grandidentata and Populus tremuloides), birch (Betula papyrifera), and cherry (Prunus pensylvatica) (???, Gates 1930, Friedman and Reich 2005, Bergen and Dronova 2007).
  • Such aspen-dominated secondary forests are sustained by background levels of disturbance, particularly fire (Gates 1930)
  • However, due to aggressive fire suppression in the US, these forests are now undergoing succession, primarily to mid/late-successional hardwood (typified by sugar maple and American beech, Fagus grandifolia) or upland conifer (typified by red and white pine) (Bergen and Dronova 2007).
  • The combination of stand-replacing disturbance ~100 years ago and suppressed disturbance since then makes this an ideal site for testing a vegetation model run from bare ground without having to prescribe disturbance.
• Guiding questions:
  • Which of these processes do we need to be able to accurately model what happens (“community succession and C cycle dynamics”) to Midwestern forests?
  • What is the cost of considering these processes, in terms of additional parametric uncertainty?
  • What are the relative contributions of parametric uncertainty (data limitation) vs. structural uncertainty (theoretical limitation?)?
• Study overview:
  • ED with factorial combination of submodels related to radiative transfer formulation (two-stream vs. multiple scatter), horizontal competition (finite canopy radius vs. complete shading), and trait plasticity (whether or not SLA and Vcmax vary with light level)
  • Ensemble analysis to look at parameter uncertainty

\[ \lambda = \frac{a_{proj}}{\mu} \]

Methods

Site description

• Secondary successional forest – bigtooth aspen (Populus grandidentata), northern red oak (Quercus rubra), red maple (Acer rubrum), paper birch (Betula papyrifera), eastern white pine (Pinus strobus)
• Average age – 95 years (2013)
• NEP in UMBS tower is 1.58 MG C ha^-1 yr^-1 (0.80-1.98) (1999-2006), but with landscape variation (REF?)
• Heavily logged in late 1800s and early 1900s, disturbed by fire until 1923 (REF)
• Present day composition typical of upper Great Lakes region
• Previously, the site of the FASET experiment (Gough et al. 2008, 2013)
• UMBS is 87% upland and 13% wetland (Bergen and Dronova 2007). We focus on the upland ecosystem, which has the following characteristics (from Bergen and Dronova (2007)):
  • 20.4% moraine, 37.8% high outwash plain, 31.3% low outwash plain, 5.7% lake-margin terrace, 3.6% ice-contact, 1.2% lowland glacial lake
  • 60.9% aspen-dominated, 16.6% northern hardwood (Acer saccharum, Acer rubrum, Fagus grandifolia, Tilia americana (basswood), Betula alleghaniensis (yellow birch), Fraxinus americana (white ash), Tsuga canadensis), 13.3% upland conifer (Pinus strobus, Pinus resinosa)
  • 1.6% northern red oak (Quercus rubra)
  • 60.9% transitioning from aspen to next succession, of which: 55.8% to northern hardwood, 34.7% to upland conifer, 3.8% to northern red oak, 2.6% to northern red oak-red maple, 3.1% to upland conifer-northern hardwood

PFT definitions

• Early successional hardwood
  • Betula papyrifera (white/paper birch)
  • Populus grandidentata (bigtooth aspen)
  • Populus tremuloides (trembling aspen)
• Mid successional hardwood
  • Quercus rubra (red oak)
  • Acer rubrum (red maple)
  • Acer pensylvaticum (striped maple)
• Late successional hardwood
  • Acer saccharum (sugar maple)
  • Fagus grandifolia (American beech)
• Pine
  • Pinus strobus (white pine)
• Omitted candidates
  • Early hardwood
    • Prunus pensylvatica (fire cherry)
  • Mid hardwood
    • Betula alleghaniensis (yellow birch)
    • Tilia americana (basswood)
    • Fraxinus americana (white ash)
  • Late hardwood
    • Ostyra virginiana (hop-hornbeam)
  • Pine
    • Pinus resinosa (red pine)

pagebreak

Colophon

This report was generated on 2019-07-23 16:38:42 using the following computational environment and dependencies:

#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.6.1 (2019-07-05)
#>  os       macOS High Sierra 10.13.6   
#>  system   x86_64, darwin17.7.0        
#>  ui       unknown                     
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2019-07-23                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package       * version    date       lib source        
#>  assertthat      0.2.1      2019-03-21 [1] CRAN (R 3.6.0)
#>  backports       1.1.4      2019-04-10 [1] CRAN (R 3.6.0)
#>  base64url       1.4        2018-05-14 [1] CRAN (R 3.6.0)
#>  callr           3.3.0      2019-07-04 [1] CRAN (R 3.6.1)
#>  cli             1.1.0      2019-03-19 [1] CRAN (R 3.6.0)
#>  codetools       0.2-16     2018-12-24 [3] CRAN (R 3.6.1)
#>  colorspace      1.4-1      2019-03-18 [1] CRAN (R 3.6.0)
#>  cowplot       * 1.0.0      2019-07-11 [1] CRAN (R 3.6.1)
#>  crayon          1.3.4      2017-09-16 [1] CRAN (R 3.6.0)
#>  data.table    * 1.12.2     2019-04-07 [1] CRAN (R 3.6.0)
#>  desc            1.2.0      2018-05-01 [1] CRAN (R 3.6.0)
#>  devtools        2.1.0      2019-07-06 [1] CRAN (R 3.6.1)
#>  digest          0.6.20     2019-07-04 [1] CRAN (R 3.6.1)
#>  dplyr         * 0.8.3      2019-07-04 [1] CRAN (R 3.6.1)
#>  drake         * 7.4.0      2019-06-07 [1] CRAN (R 3.6.0)
#>  evaluate        0.14       2019-05-28 [1] CRAN (R 3.6.0)
#>  forcats       * 0.4.0      2019-02-17 [1] CRAN (R 3.6.0)
#>  fortebaseline * 0.0.0.9000 2019-07-23 [1] local         
#>  fs            * 1.3.1      2019-05-06 [1] CRAN (R 3.6.0)
#>  fst           * 0.9.0      2019-04-09 [1] CRAN (R 3.6.0)
#>  furrr         * 0.1.0      2018-05-16 [1] CRAN (R 3.6.0)
#>  future        * 1.14.0     2019-07-02 [1] CRAN (R 3.6.1)
#>  future.callr  * 0.4.0      2019-01-07 [1] CRAN (R 3.6.0)
#>  ggplot2       * 3.2.0      2019-06-16 [1] CRAN (R 3.6.0)
#>  ggrepel       * 0.8.1      2019-05-07 [1] CRAN (R 3.6.1)
#>  globals         0.12.4     2018-10-11 [1] CRAN (R 3.6.0)
#>  glue            1.3.1      2019-03-12 [1] CRAN (R 3.6.0)
#>  gtable          0.3.0      2019-03-25 [1] CRAN (R 3.6.0)
#>  here          * 0.1        2017-05-28 [1] CRAN (R 3.6.0)
#>  highr           0.8        2019-03-20 [1] CRAN (R 3.6.0)
#>  hms             0.5.0      2019-07-09 [1] CRAN (R 3.6.1)
#>  htmltools       0.3.6      2017-04-28 [1] CRAN (R 3.6.0)
#>  igraph          1.2.4.1    2019-04-22 [1] CRAN (R 3.6.0)
#>  knitr           1.23       2019-05-18 [1] CRAN (R 3.6.0)
#>  labeling        0.3        2014-08-23 [1] CRAN (R 3.6.0)
#>  lazyeval        0.2.2      2019-03-15 [1] CRAN (R 3.6.0)
#>  listenv         0.7.0      2018-01-21 [1] CRAN (R 3.6.0)
#>  lubridate     * 1.7.4      2018-04-11 [1] CRAN (R 3.6.0)
#>  magrittr        1.5        2014-11-22 [1] CRAN (R 3.6.0)
#>  memoise         1.1.0      2017-04-21 [1] CRAN (R 3.6.0)
#>  munsell         0.5.0      2018-06-12 [1] CRAN (R 3.6.0)
#>  pillar          1.4.2      2019-06-29 [1] CRAN (R 3.6.0)
#>  pkgbuild        1.0.3      2019-03-20 [1] CRAN (R 3.6.0)
#>  pkgconfig       2.0.2      2018-08-16 [1] CRAN (R 3.6.0)
#>  pkgload         1.0.2      2018-10-29 [1] CRAN (R 3.6.0)
#>  plyr            1.8.4      2016-06-08 [1] CRAN (R 3.6.0)
#>  prettyunits     1.0.2      2015-07-13 [1] CRAN (R 3.6.0)
#>  processx        3.4.0      2019-07-03 [1] CRAN (R 3.6.1)
#>  ps              1.3.0      2018-12-21 [1] CRAN (R 3.6.0)
#>  purrr         * 0.3.2      2019-03-15 [1] CRAN (R 3.6.0)
#>  R6              2.4.0      2019-02-14 [1] CRAN (R 3.6.0)
#>  Rcpp            1.0.1      2019-03-17 [1] CRAN (R 3.6.0)
#>  readr         * 1.3.1      2018-12-21 [1] CRAN (R 3.6.0)
#>  remotes         2.1.0      2019-06-24 [1] CRAN (R 3.6.1)
#>  reshape2        1.4.3      2017-12-11 [1] CRAN (R 3.6.0)
#>  rlang           0.4.0      2019-06-25 [1] CRAN (R 3.6.0)
#>  rmarkdown       1.14       2019-07-12 [1] CRAN (R 3.6.1)
#>  rprojroot       1.3-2      2018-01-03 [1] CRAN (R 3.6.0)
#>  scales          1.0.0      2018-08-09 [1] CRAN (R 3.6.0)
#>  sessioninfo     1.1.1      2018-11-05 [1] CRAN (R 3.6.0)
#>  storr           1.2.1      2018-10-18 [1] CRAN (R 3.6.0)
#>  stringi         1.4.3      2019-03-12 [1] CRAN (R 3.6.0)
#>  stringr         1.4.0      2019-02-10 [1] CRAN (R 3.6.0)
#>  testthat        2.1.1      2019-04-23 [1] CRAN (R 3.6.0)
#>  tibble          2.1.3      2019-06-06 [1] CRAN (R 3.6.0)
#>  tidyr         * 0.8.3      2019-03-01 [1] CRAN (R 3.6.0)
#>  tidyselect      0.2.5      2018-10-11 [1] CRAN (R 3.6.0)
#>  usethis         1.5.1      2019-07-04 [1] CRAN (R 3.6.1)
#>  vctrs           0.2.0      2019-07-05 [1] CRAN (R 3.6.1)
#>  withr           2.1.2      2018-03-15 [1] CRAN (R 3.6.0)
#>  xfun            0.8        2019-06-25 [1] CRAN (R 3.6.1)
#>  yaml            2.2.0      2018-07-25 [1] CRAN (R 3.6.0)
#>  zeallot         0.1.0      2018-01-28 [1] CRAN (R 3.6.0)
#> 
#> [1] /Users/shik544/R/library
#> [2] /usr/local/lib/R/3.6/site-library
#> [3] /usr/local/Cellar/r/3.6.1_1/lib/R/library

The current Git commit details are:

#> Local:    master /Users/shik544/local_projects/fortebaseline
#> Remote:   master @ origin (https://github.com/ashiklom/fortebaseline)
#> Head:     [60dfe29] 2019-07-23: Update sensitivity analysis code and text

View this project on GitHub in repository ashiklom/fortebaseline (this manuscript file is in analysis/paper/paper.Rmd) and in Open Science Framework (OSF) at https://osf.io/dznuf/ .